Poravnanja - 1. deo

Čest zadatak bioinformatičara je pronalaženje poravnanja dve ili više sekvenci. Poravnanje podrazumeva pronalaženje najboljeg poklapanja karaktera sekvenci sa ciljem pronalaženja regiona velike sličnosti. Uzmimo, na primer, aminokiselinsku FASTA sekvencu histonskog proteina H1 kod čoveka:

>sp|P07305-2|H10_HUMAN Isoform 2 of Histone H1.0 OS=Homo sapiens OX=9606 GN=H1-0
MTKKSTDHPKYSDMIVAAIQAEKNRAGSSRQSIQKYIKSHYKVGENADSQIKLSIKRLVT
TGVLKQTKGVGASGSFRLAKSDEPKKSVAFKKTKKEIKKVATPKKASKPKKAASKAPTKK
PKATPVKKAKKKLAATPKKAKKPKTVKAKPVKASKPKKAKPVKPKAKSSAKRAGKKK

i uporedimo je sa sekvencom istog proteina kod miša:

>sp|P10922|H10_MOUSE Histone H1.0 OS=Mus musculus OX=10090 GN=H1-0 PE=2 SV=4
MTENSTSAPAAKPKRAKASKKSTDHPKYSDMIVAAIQAEKNRAGSSRQSIQKYIKSHYKV
GENADSQIKLSIKRLVTTGVLKQTKGVGASGSFRLAKGDEPKRSVAFKKTKKEVKKVATP
KKAAKPKKAASKAPSKKPKATPVKKAKKKPAATPKKAKKPKVVKVKPVKASKPKKAKTVK
PKAKSSAKRASKKK

Sekvence ova dva proteina nisu potpuno identične, ali se primećuje velika sličnost. Poravnanje ove dve sekvence, koje bi prepoznalo regione velike sličnosti moglo bi izgledati ovako:

H10_HUMAN          1 -----------------MTKKSTDHPKYSDMIVAAIQAEKNRAGSSRQSI     33
                                      .:|||||||||||||||||||||||||||||||
H10_MOUSE          1 MTENSTSAPAAKPKRAKASKKSTDHPKYSDMIVAAIQAEKNRAGSSRQSI     50

H10_HUMAN         34 QKYIKSHYKVGENADSQIKLSIKRLVTTGVLKQTKGVGASGSFRLAKSDE     83
                     |||||||||||||||||||||||||||||||||||||||||||||||.||
H10_MOUSE         51 QKYIKSHYKVGENADSQIKLSIKRLVTTGVLKQTKGVGASGSFRLAKGDE    100

H10_HUMAN         84 PKKSVAFKKTKKEIKKVATPKKASKPKKAASKAPTKKPKATPVKKAKKKL    133
                     ||:||||||||||:|||||||||:||||||||||:||||||||||||||.
H10_MOUSE        101 PKRSVAFKKTKKEVKKVATPKKAAKPKKAASKAPSKKPKATPVKKAKKKP    150

H10_HUMAN        134 AATPKKAKKPKTVKAKPVKASKPKKAKPVKPKAKSSAKRAGKKK    177
                     |||||||||||.||.||||||||||||.||||||||||||.|||
H10_MOUSE        151 AATPKKAKKPKVVKVKPVKASKPKKAKTVKPKAKSSAKRASKKK    194

Pri čemu uspravne crte | označavaju uspešno uparivanje karaktera, horizontalne crtice - uparivanje sa prazninom (insercije/delecije) dok tačkice :. označavaju manju ili veću grešku prilikom uparivanja različitih karaktera.

Iz navedenih rezultata poravnanja može se zaključiti da je H1 protein ostao većinski nepromenjen tokom evolucije ove dve vrsta organizama.

Kao što se može primetiti, na algoritmima poravnanja je zadatak poravnati slične, ne nužno identične, sekvence. Iz ove nedovoljno precizne definicije zadatka mogu se na različite načine konstruisati algoritmi sa različitim ciljevima poravnanja.

Većina ovih algoritama su predstavnici algoritama dinamičkog programiranja pa će iz tog razloga prvo biti definisan jednostavniji problem, kako bi se bolje shvatila problematika i način rešavanja problema.

1. Manhattan Tourist Problem

Kao primer algoritma dinamičkog programiranja, biće predstavljen algoritam koji rešava Manhattan Tourist Problem. Zadatak je pomoći turisti da obiđe ulice Menhetna krećući se od levog ugla mape ka donjem desnom uglu. Kako su na određenim ulicama atrakcije (čija je zanimljivost definisana brojem), cilj turiste je da poseti što više zanimljivih atrakcija tako da na kraju obilaska ima najveći mogući ukupan zbir zanimljivosti posećenih lokacija. Atrakcija se smatra posećenom ako je turista prošao ulicom u kojoj se atrakcija nalazi.

Rešenje

Svaki problem dinamičkog progamiranja, kao i njegovo rešenje, mogu se predstaviti direktnim acikličnim grafom (DAG). Svaki početak ulice i raskrsnica mogu se predstaviti čvorom grafa sa oznakom $i \cdot m + j$, gde $i$ i $j$ predstavljaju raskrsnice $i$-te i $j$-te ulice. Vrednosti $i=0$ i $j = 0$ označavaju početke odgovarajućih ulica. Zanimljivosti (atrakcija koje se nalaze u ulici) predstavljene su granama grafa, između čvorova, sa težinama koje odgovaraju zanimljivosti atrakcija u delu ulice između navedenih raskrsnica. Čvorovi grafa se mogu predstaviti matricom dimenzija $(n+1)\textrm{ }x\textrm{ }(m+1)$ dok se grane mogu predstaviti uređenim trojkama $(u,v,w)$, gde su $u$ i $v$ čvorovi između kojih se nalazi grana a $w$ težina grane.

In [1]:
'''
Mapa grada, [x] su čvorovi grafa a -y- su težine grana
[0] - 1 - [1] - 1 - [2]
|          |         |
1          2         1
|          |         |
[3] - 1 - [4] - 1 - [5]
|          |         |
1          2         1
|          |         |
[6] - 1 - [7] - 2 - [8]

'''

n = 3
m = 3
edges = [(0,1,1),(1,2,1),(0,3,1),(1,4,2),(2,5,1),(3,4,1),(4,5,1),(3,6,1),(4,7,2),(5,8,1),(6,7,1),(7,8,2)]

Kako nije efikasno prolaziti kroz listu grana svaki put kada je potrebno proveriti težinu, listu možemo pretvoriti u odgovarajuću mapu $[u]\rightarrow[v]\rightarrow[w]$

In [2]:
def edges_to_map(edges):
    edge_map = {}
    for (u, v, w) in edges:
        if u not in edge_map:
            edge_map[u] = {}
        edge_map[u][v] = w

    return edge_map
In [3]:
edges_to_map(edges)
Out[3]:
{0: {1: 1, 3: 1},
 1: {2: 1, 4: 2},
 2: {5: 1},
 3: {4: 1, 6: 1},
 4: {5: 1, 7: 2},
 5: {8: 1},
 6: {7: 1},
 7: {8: 2}}

Turista se može kretati samo u pravcima dole i desno. Za svaku poziciju u grafu može se izračunati cena optimalnog puta (sa najviše posećenih atrakcija) rekurentnom formulom: $$Score_{i,j} = max\begin{cases} Score_{i-1,j} + edge_{i-1,j}\\ Score_{i,j-1} + edge_{i,j-1} \end{cases}$$

Drugim rečima, do polja ${i,j}$ može se doći odozgo, sa polja ${i-1,j}$ - čime se na skor sakupljen do te pozicije dodaje vrednost grane od $i-1,j$ do $i,j$, ili sa leve strane, sa polja ${i,j - 1}$ čime se na skor sakupljen do te pozicije dodaje vrednost grane od $i,j - 1$ do $i,j$. Pri takvom izboru, za novi skor bira se maksimum ponuđenih mogućnosti. Do polja koja se nalaze na gornjoj ivici ($i = 0$) može se doći samo kretanjem sa leve strane pa je cena dolaska sa gornje strane uvek jednaka nuli. Analogno važi i za levu ivicu ($j + 0$). Iz tog razloga, matricu $Score$ možemo proširiti dodatnim redom i kolonom sa gornje i leve strane, čije će vrednosti biti jednake nuli a obilazak će počinjati od pozicije ($i = 1, $j = 1). Mapiranje koordinata u čvorove grafa će sada biti $(i-1) * m + (j-1)$. Vrednost na poziciji $(n + 1, m + 1)$ će biti vrednost maksimalnog skora atrakcija pri obilasku dok će odabir puta, pri svakom izboru na putu od početka do kraja obilaska, odgovarati koracima optimalnom obilaska grada.

In [4]:
def manhattan_tourist_problem(m, n, edges):
    S = [[0 for j in range(m + 1)] for i in range(n + 1)]
    backtrack = [[0 for j in range(m + 1)] for i in range(n + 1)]
    edge_map = edges_to_map(edges)
    
    # Popunjavanje tabele dinamickog programiranja
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            u_up = (i - 2) * m + (j - 1)
            u_left = (i - 1) * m + (j - 2)
            v = (i - 1) * m + (j - 1)
            
                #          u_up 
                #           |
                # u_left -- v
                # 
             
            top_edge = 0
            if u_up in edge_map:
                if v in edge_map[u_up]:
                    top_edge = edge_map[u_up][v]
                
            left_edge = 0
            if u_left in edge_map:
                if v in edge_map[u_left]:
                    left_edge = edge_map[u_left][v]
                
            from_top = S[i - 1][j] + top_edge
            from_left = S[i][j - 1] + left_edge
            
            S[i][j] = max(from_top, from_left)
            
            # Provera koji put je odabran do tekuce pozicije
            if S[i][j] == from_top:
                backtrack[i][j] = 1
            else:
                backtrack[i][j] = -1
                
            
    i = n
    j = m
    
    path = [(i-1) * m + (j - 1)]
    
    while i != 1 or j != 1:
        if backtrack[i][j] == 1:
            path = [(i - 2) * m + (j-1)] + path
            i -= 1
        else:
            path = [(i-1) * m + (j - 2)] + path
            j -= 1
            
        
        
    return S[n][m], path
In [5]:
(score, path) = manhattan_tourist_problem(m, n, edges)
print(f'Ukupan skor: {score}')
print(f'Putanja: {path}')
Ukupan skor: 7
Putanja: [0, 1, 4, 7, 8]

2. Najduža zajednička podsekvenca

Za početak, treba napraviti razliku između podniske i podsekvence. Poniska podrazumeva zajednički niz uzastopnih karaktera koji se nalazi u okviru niske dok podsekvenca podrazumeva niz karaktera koji ne moraju nužno biti uzastopni, ali im je poredak isti kao u originalnoj niski.

Primer

Niska: ABCDEFG

Podniske: ABC, CDEF, FG, ...

Podsekvence: ACG, ABFG, CE, ...

Najduža zajednička podsekvenca između dve niske je najduža podsekvenca koja se nalazi i u jednoj u u drugoj niski. Problem se rešava dinamičkim programiranjem. Pretpostavimo da je poznata dužina podsekvence do prefiksa $i-1$ u prvoj niski i $j$ u drugoj niski i prefiksa $i$ u prvoj i $j-1$ u drugoj niski. Najduža podsekvenca do prefiksa $i$ i $j$ produžuje duži od prethodna dva prefiksa za jedan, ako su karakteri $niska1_i$ i $niska2_j$ jednaki ili ostaje skor dužeg prefiksa ukoliko se karakteri razlikuju. Vrednosti najdužih podsekvenci svih prefiksa obe niske mogu se predstaviti tabelom, tako da se na poziciji $i,j$ nalazi dužina najduže podsekvence prefiksa $i$ karaktera prve niske i $j$ karaktera druge niske. Tabela će biti dimenzija $(n+1)\textrm{ }x\textrm{ }(m+1)$, gde su $n$ i $m$ dužine niski, pri čemu nulti red i nulta kolona odgovaraju praznim prefiksima (dužine 0) odgovarajućih niski. Za bazu rekurzije može se uzeti činjenica da poravnanje praznog prefiksa sa bilo čime daje najdužu podsekvencu dužine 0. Iz toga sledi da su vrednosti u redu $i = 0$ i koloni $j = 0$ jednake nuli. Vrednost u polju $n,m$ popiunjene tabele odgovara odgovoru na polazni problem - dužini najduže zajedniče podsekvence.

$$Score_{i,j} = max\begin{cases} Score_{i-1,j-1} + 1,& niska_i = niska_j\\ Score_{i,j-1},& niska_i \neq niska_j\\ Score_{i,j-1},& niska_i = niska_j\\ \end{cases}$$
In [6]:
def lcs_backtrack(v, w):
    n = len(v) + 1
    m = len(w) + 1
    
    s = [[0 for j in range(m)] for i in range(n)]
    backtrack = [[None for j in range(m)] for i in range(n)]
    
    for i in range(1, n):
        backtrack[i][0] = (i - 1, 0)
        
    for j in range(1, m):
        backtrack[0][j] = (0, j - 1)
    
    for i in range(1, n):
        for j in range(1, m):
            from_top = s[i - 1][j] + 0
            from_left = s[i][j - 1] + 0
            from_diagonal = s[i - 1][j - 1] + int(v[i - 1] == w[j - 1])
            
            s[i][j] = max(from_top, from_left, from_diagonal)
            
            # Opciono, pamćenje pozicija sa kojih je napravljen korak ka tekućoj
            if s[i][j] == from_top:
                backtrack[i][j] = (i - 1, j)
            elif s[i][j] == from_left:
                backtrack[i][j] = (i, j - 1)
            else:
                backtrack[i][j] = (i - 1, j - 1)
                
    i = n - 1
    j = m - 1
    
    lcs = ""
    
    # Opciono, rekonstrukcija najduže zajedničke podsekvence
    while backtrack[i][j] != None:
        if backtrack[i][j] == (i - 1, j - 1):
            lcs = v[i - 1] + lcs
        (i, j) = backtrack[i][j]
           
    length = s[n - 1][m - 1]
    
    return length, lcs
In [7]:
v = 'AGTCGTGATCGTTGTA'
w = 'GTATGAA'

length, lcs = lcs_backtrack(v, w)

print('Dužina najduže zajedničke podsekvence: {}'.format(length))
print('Najduža zajednička podsekvenca: {}'.format(lcs))
Dužina najduže zajedničke podsekvence: 6
Najduža zajednička podsekvenca: GTATGA

3. Edit rastojanje

Edit rastojanje između dve niske predstavlja minimalni broj edit operacija (dodavanje karaktera, brisanje karaktera, zamena karaktera) potrebnih da se jedna niska transformiše u drugu.

Primer:

Niska 1: 'ABC' Niska 2: 'ABD'

Edit rastojanje između navedenih niski je 1 (jedna zamena karaktera)

Niska 1: 'AB' Niska 2: 'ABC'

Edit rastojanje između navedenih niski je takođe 1 (jedno dodavanje karaktera)

Niska 1: 'ABC' Niska 2: 'AD'

Edit rastojanje između navedenih je takođe 2 (jedno brisanje karaktera i jedna zamena)

Ovaj algoritam je takođe predstavnik algoritama dinamičkog programiranja i može se koristiti za potrebe poravnanja bioloških sekvenci. Konstrukcija algoritma je ponovo rekurzivna, Neka je poznato edit rastojanje prefiksa prve niske do karaktera na poziciji $i-1$ i karaktera na poziciji $j-1$ druge niske. Edit rastojanje prefiksa $i$ i $j$ može se dobiti na tri načina:

  • karakteri na pozicijama niska1[i] i niska2[j] su jednaki $\rightarrow$ Edit rastojanje je jednako edit rastojanju prefiksa $i-1$ i $j-1$ (cena je 0 jer nije potrebno primeniti ni jednu dodatnu edit operaciju) ili se karakteri na pozicijama niska1[i] i niska2[j] se razlikuju $\rightarrow$ Potrebno je izvršiti još jednu edit operaciju kako bi se karakteri zamenili

  • niska1[i] se poravnava sa niska2[j-1] $\rightarrow$ izvršeno je brisanje karaktera iz niske 1 (još jedna edit operacija)

  • niska1[i-1] se poravnava sa niska2[j] $\rightarrow$ izvršeno je brisanje karaktera iz niske 2 (još jedna edit operacija)

Za poravnanje prazne sekvence sa sekvencom duzine x potrebno je x edit operacija. Upravo će to biti baza rekurzije. Potrebno je konstruisati tabelu dinamičkog programiranja dimenzije $(n+1)\textrm{ }x\textrm{ }(m+1)$, gde su $n$ i $m$ dužine niski. Dodatak $+1$ podrazumeva da nulti red i nulta kolona predstavljaju poravnanja sa praznom niskom. U svakom koraku, cilj je odabrati ponuđeni put koji vodi ka minimalnom edit rastojanju.

Red $i = 0$ i kolona $j = 0$ odgovaraju poravnanjima sa praznom niskom, prema tome vrednost u tim poljima će biti $S_{0,j} = j$ i $S_{i,0} = i$.

Ostatak tabele popunjava se na osnovu rekurentne jednačine:

$$Score_{i,j} = min\begin{cases} Score_{i-1,j-1} + 0,& niska_i = niska_j\\ Score_{i-1,j-1} + 1,& niska_i \neq niska_j\\ Score_{i-1,j} + 1\\ Score_{i,j-1} + 1 \end{cases}$$
In [8]:
def edit_distance(v, w):
    n = len(v) + 1
    m = len(w) + 1
    
    s = [[0 for j in range(m)] for i in range(n)]
    backtrack = [[None for j in range(m)] for i in range(n)]
    
    for i in range(1, n):
        s[i][0] = i
        backtrack[i][0] = (i - 1, 0)
        
    for j in range(1, m):
        s[0][j] = j
        backtrack[0][j] = (0, j - 1)
        
    for i in range(1, n):
        for j in range(1, m):

            from_top = s[i - 1][j] + 1
            from_left = s[i][j - 1] + 1
            from_diagonal = s[i - 1][j - 1] + int(v[i - 1] != w[j - 1])
            
            s[i][j] = min(from_top, from_left, from_diagonal)
            
            # Opciono, rekonstrukcija poravnanja
            if s[i][j] == from_top:
                backtrack[i][j] = (i - 1, j)
            elif s[i][j] == from_left:
                backtrack[i][j] = (i, j - 1)
            else:
                backtrack[i][j] = (i - 1, j - 1)
                
    i = n - 1
    j = m - 1
    
    v_align = ''
    w_align = ''
    
    # Opciono, rekonstrukcija poravnanja
    while backtrack[i][j] != None:
        if backtrack[i][j] == (i - 1, j):
            v_align = v[i - 1] + v_align 
            w_align = '-' + w_align
        elif backtrack[i][j] == (i, j - 1):
            v_align = '-' + v_align 
            w_align = w[j - 1] + w_align
        else:
            v_align = v[i - 1] + v_align
            w_align = w[j - 1] + w_align
            
        (i, j) = backtrack[i][j]
    
    distance = s[n - 1][m - 1]
    
    return distance, v_align, w_align
In [9]:
v = 'ACGTGACCTGGA'
w = 'GTGT'

distance, v_align, w_align = edit_distance(v, w)

print('Edit distance: {}'.format(distance))
print(v_align)
print(w_align)
Edit distance: 8
ACGTGACCTGGA
--GTG---T---

4. Globalno poravnanje

Ponekad edit rastojanje nije dovoljno da bi objasnilo određena poravnanja već je potrebno dodeliti određeni skor, odnosno kaznu, za odgovarajuća poravnanja karaktera. Osnova za tako nešto može se naći u biološkim osobinama proteinskih sekvenci. Aminokiseline imaju odgovarajuća svojstva koja utiču na način na koji se protein uvija u svoju 3D strukturu. Ako bi se između dve sekvence poravnale aminokiseline sa potpuno drugačijim svojstvima, kazna takvog poravnanja bila bi veća nego u slučaju poravnanja dve aminokiseline sa sličnim svojstvima. Takođe kazna pojave praznine u poravnanju ne mora biti jednaka kazni loše uparenih aminokiselina. Iz tih razloga, pojavljuje se pojam globalnog poravnanja sa različitim kaznama za loša poravnanja i praznine i pozitivnim skorovima za dobra poravnanja. Za razliku od edit rastojanja, koje računa minimum rastojanja prethodnih prefiksa, globalno poravnanje (Needleman-Wunsch algoritam) računa maksimum ukupnog skora. Ukupan skor poravnanja dve niske nalazi se u polju $n+1,m+1$ tabele.

$$Score_{0,j} = j * GapScore$$$$Score_{i,0} = i * GapScore$$$$Score_{i,j} = max\begin{cases} Score_{i-1,j-1} + GoodMatchScore,& niska_i = niska_j\\ Score_{i-1,j-1} + MissMatchScore,& niska_i \neq niska_j\\ Score_{i-1,j} + GapScore\\ Score_{i,j-1} + GapScore \end{cases}$$
In [10]:
GAP_PENALTY = -2
MISSMATCH_PENALTY = 0
MATCH_SCORE = 1

def needleman_wunsch(v, w):
    n = len(v) + 1
    m = len(w) + 1
    
    s = [[0 for j in range(m)] for i in range(n)]
    backtrack = [[None for j in range(m)] for i in range(n)]
    
    for i in range(1, n):
        s[i][0] = s[i - 1][0] + GAP_PENALTY
        backtrack[i][0] = (i - 1, 0)
        
    for j in range(1, m):
        s[0][j] = s[0][j - 1] + GAP_PENALTY
        backtrack[0][j] = (0, j - 1)
        
    for i in range(1, n):
        for j in range(1, m):

            from_top = s[i - 1][j] + GAP_PENALTY
            from_left = s[i][j - 1] + GAP_PENALTY
            
            if v[i - 1] == w[j - 1]:
                from_diagonal = s[i - 1][j - 1] + MATCH_SCORE
            else:
                from_diagonal = s[i - 1][j - 1] + MISSMATCH_PENALTY
                
            s[i][j] = max(from_top, from_left, from_diagonal)
            
            # Opciono, rekonstrukcija poravnanja
            if s[i][j] == from_top:
                backtrack[i][j] = (i - 1, j)
            elif s[i][j] == from_left:
                backtrack[i][j] = (i, j - 1)
            else:
                backtrack[i][j] = (i - 1, j - 1)
                
    i = n - 1
    j = m - 1
    
    v_align = ''
    w_align = ''
    
    # Opciono, rekonstrukcija poravnanja
    while backtrack[i][j] != None:
        if backtrack[i][j] == (i - 1, j):
            v_align = v[i - 1] + v_align 
            w_align = '-' + w_align
        elif backtrack[i][j] == (i, j - 1):
            v_align = '-' + v_align 
            w_align = w[j - 1] + w_align
        else:
            v_align = v[i - 1] + v_align
            w_align = w[j - 1] + w_align
            
        (i, j) = backtrack[i][j]
    
    score = s[n - 1][m - 1]
    
    return score, v_align, w_align
In [11]:
v = 'ACGTGACCTGGA'
w = 'GTGT'

score, v_align, w_align = needleman_wunsch(v, w)

print('Skor globalnog poravnanja: {}'.format(score))
print(v_align)
print(w_align)
Skor globalnog poravnanja: -12
ACGTGACCTGGA
--GTG---T---

5. Lokalno poravnanje

Globalno poravnanje pokušava da poravna cele sekvence jednu sa drugom. U nekim situacijama potrebno je samo pronaći lokalna poravnanja sa najvećim skorom (najbolje očuvane regione). Modifikacija globalnog poravnanja u lokalno je jednostavna. Baza rekurzije podrazumeva da vrednosti reda $i=0$ i $j=0$ budu 0, kako se ne bi kažnjavale početne praznine dok se za rezultat lokalnog poravnanja ne uzima vrednost iz tabele na poziciji $n+1,m+1$ već maksimalna vrednost u tabeli.

$$Score_{0,j} = 0$$$$Score_{i,0} = 0$$$$Score_{i,j} = max\begin{cases} Score_{i-1,j-1} + GoodMatchScore,& niska_i = niska_j\\ Score_{i-1,j-1} + MissMatchScore,& niska_i \neq niska_j\\ Score_{i-1,j} + GapScore\\ Score_{i,j-1} + GapScore \end{cases}$$$$LocalAlignmentScore = max(S)$$
In [12]:
GAP_PENALTY = -1
MISSMATCH_PENALTY = -2
MATCH_SCORE = 1

def smith_waterman(v, w):
    n = len(v) + 1
    m = len(w) + 1
    
    s = [[0 for j in range(m)] for i in range(n)]
    backtrack = [[None for j in range(m)] for i in range(n)]
    
    for i in range(1, n):
        backtrack[i][0] = (i - 1, 0)
        
    for j in range(1, m):
        backtrack[0][j] = (0, j - 1)
        
    for i in range(1, n):
        for j in range(1, m):

            from_top = s[i - 1][j] + GAP_PENALTY
            from_left = s[i][j - 1] + GAP_PENALTY
            
            if v[i - 1] == w[j - 1]:
                from_diagonal = s[i - 1][j - 1] + MATCH_SCORE
            else:
                from_diagonal = s[i - 1][j - 1] + MISSMATCH_PENALTY
                
            s[i][j] = max(from_top, from_left, from_diagonal, 0)
            
            # Opciono, rekonstrukcija poravnanja
            if s[i][j] == from_top:
                backtrack[i][j] = (i - 1, j)
            elif s[i][j] == from_left:
                backtrack[i][j] = (i, j - 1)
            else:
                backtrack[i][j] = (i - 1, j - 1)
            
    max_value = float('-inf')        
    
    for i in range(n):
        for j in range(m):
            if s[i][j] > max_value:
                max_value = s[i][j]
                max_i = i
                max_j = j
            
    i = max_i
    j = max_j
    
    v_align = ''
    w_align = ''
    
    # Opciono, rekonstrukcija poravnanja
    while s[i][j] != 0:
        if backtrack[i][j] == (i - 1, j):
            v_align = v[i - 1] + v_align 
            w_align = '-' + w_align
        elif backtrack[i][j] == (i, j - 1):
            v_align = '-' + v_align 
            w_align = w[j - 1] + w_align
        else:
            v_align = v[i - 1] + v_align
            w_align = w[j - 1] + w_align
            
        (i, j) = backtrack[i][j]
    
    v_align_len = len(v_align)
    w_align_len = len(w_align)
    
    if m < n:
        w_align = (i * '-') + w + (n - max_i - 2) * '-'
        v_align = v
    else:
        v_align = (j * '-') + v + (m - max_j) * '-'
        w_align = w
    
    score = s[n - 1][m - 1]
    
    return score, v_align, w_align
In [13]:
v = 'ACGTGACCTGGA'
w = 'GTGGT'

score, v_align, w_align = smith_waterman(v, w)

print('Ukupan skor lokalnog poravnanja: {}'.format(score))
print(v_align)
print(w_align)
Ukupan skor lokalnog poravnanja: 1
ACGTGACCTGGA
--GTGGT------